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Abstract. We analyze mean- field models of synchronization of phase oscillators with 
singular couplings and subject to external random forces. They are related to the 
Kuramoto-Sakaguchi model. Their probability densities satisfy local partial differential 
equations similar to the Porous Medium, Burgers and extended Burgers equations 
depending on the degree of singularity of the coupling. We show that Porous Medium 
oscillators (the most singularly coupled) do not synchronize and that (transient) 
synchronization is possible only at zero temperature for Burgers oscillators. The 
extended Burgers oscillators have a nonlocal coupling first introduced by Daido and 
they may synchronize at any temperature. Exact expressions for their synchronized 
phases and for Daido's order function are given in terms of elliptic functions. 



1. Introduction 

Collective synchronization of large populations of nonlinearly coupled phase oscillators 
has been intensely studied since Winfree realized its importance for biological systems [Q] 
and Kuramoto gave mathematical form to these ideas in a simple model 0, ^. 
Motivation for studying the Kuramoto model can be found in the broad variety of 
physical, chemical or biological phenomena which can be modelled within its framework 
(see D, and references therein). 

The problem we want to consider is the dynamics of a system of nonlinearly, globally 
coupled phase oscillators with random frequencies uji [the probability for an oscillator to 
have frequency uj G {u, u+du) isp^u) du], subject to external (independent, identically- 
distributed) white noise sources rji (of strength \pXT\. 

^ = ^.-§E/(0.-0.) + ^.W, (1) 

i = \, . . . ,N . Here 0j(t) denotes the zth oscillator phase, if > represents the coupling 
strength, and / is a generic real function of periodicity In. In Fourier space, the latter 
can be decomposed as follows, 

oo 

f{^)= E anexp(ma;) (2) 

?i=— oo 

with a„ = al„. 

Important work on the Kuramoto model with a general coupling function was 



carried out by Daido P, [T^] and Crawford [|1T] among others. Daido introduced 
the concept of order function to understand synchronization at zero temperature, 
T = 0. Essentially, the order function is an average of the oscillator drift velocity 
on a rotating frame [^, |T^. Assuming a reasonable shape for the order function and 
the one-oscillator probability density (supposedly stationary in the rotating frame), 
Daido derived a functional equation for the order function. He then found solutions 
by direct numerical simulations and by bifurcation theory for stationary or rotating- 
wave probability densities. However, his theory is not sufficiently general to cover other 
possible probability densities (e.g., standing waves O, llB]) nor to predict their stability 
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properties. Crawford considered the problem of constructing one-oscillator probability 
densities which bifurcate from the incoherent density l/(27r) (for which an oscillator has 
the same probability to be at any given value of the angle 0) |]lT|. As a consequence of 
these works, different conjectures were proposed, both concerning scaling of bifurcating 
solution branches and general statements such as "adding noise could actually yield an 
enhancement in the level of synchronization" [^. One difficulty with these previous 
works is that explicit calculations were performed mostly near bifurcation points. One 
exception is Daido's coupling function f{x) = signx for which the calculation of the 
order function can be explicitly carried out. 

In this paper, we shall introduce classes of coupling functions leading to drift 
terms which are local functions of the probability density. We restrict ourselves to 
oscillators without frequency disorder, i.e., with frequency distribution p{uj) = S{uj). 
The corresponding governing equations for the probability density (in the limit of 
infinitely many oscillators) are systems of nonlinear partial differential equations which 
are analyzed. Then definite results can be proved about synchronization which do not 
depend on bifurcation calculations, order function theory or numerical simulations as 
much as the work of previous authors. Together with singular perturbation calculations, 
our results can be used to give approximations to the one-oscillator probability density 
of models with disorder, in the limiting case of high-frequencies |T^. The examples 
of singular coupling functions considered in this paper are the hard- needles (+ sign) 
and stick- needles (- sign) couplings, for which / = ±6', the Burgers coupling, f = S, 
and the Daido coupling (extended Burgers model), /(x) = sign(x). In all cases, these 
/'s are extended periodically outside the interval — vr < x < tc. We shall show 
that less singular couplings give rise to models for which the oscillators are easier to 
synchronize: (i) the hard-needles oscillators do not synchronize at any temperature 
(the stick-needles coupling leads to ill-posed problems unless the temperature is large 
enough), (ii) the Burgers oscillators may synchronize only at zero temperature; even 
then, their probability density algebraically decays to incoherence as t~^ for large times, 
(iii) the oscillators with Daido coupling may synchronize at any temperature, and the 
synchronized phases are described by explicitly-known order functions and probability 
densities. A first paper in the direction of these results was published by the authors in 

Hi- 

The rest of the paper is structured as follows. Section |^ contains results which hold 
for general models with disorder in the natural frequencies. These include a derivation 
of the governing equations for the one-oscillator probability density, its moments and 
the moment-generating function. We establish a relationship between these objects and 
Daido's order function showing that the latter is proportional to the average oscillator 
drift in a rotating frame. Lastly, we give the leading-order form of the probability density 
in the limit of high frequencies. The density is a superposition of probability densities 
with zero natural frequency in rotating frames and with smaller coupling constants. In 
Section |^ we extend Daido's approach to the case of nonzero temperature, and find 
the form of the stationary solutions and a general Liapunov functional for the case 
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of odd coupling functions. In section ^ we introduce the family of singular couplings 
to be studied. Section ^ contains our analysis of the porous Medium models which 
appear for the hard-needles and stick-needles couplings. We show that the hard-needles 
model is well-posed, its solution exists globally, and also that sharp time decay estimates 
towards its equilibrium can be obtained locally for any initial data or globally for small 
initial data evolving exponentially fast towards incoherence. The last result is proven 
by different methods according to whether we allow the temperature to be zero. On the 
other hand, the stick-needles model is ill-posed for low enough temperature. Section ^ 
is devoted to study a model with Burgers coupling. We adapt well-known results to the 
case of periodic boundary conditions to show that the probability density of the Burgers 
oscillators may tend to a state different from incoherence only if the temperature is zero. 
In section ^ we analyze the extended Burgers model resulting from the Daido coupling 
/(x) = sign(x). We first show that finding the probability density is equivalent to solving 
two coupled nonlinear parabolic equations for the drift v{x,t) and a certain functional 
of the probability density, a{x,t). Then we find families of stationary solutions in 
terms of elliptic functions. These solutions bifurcate supercritically from incoherence 
at couplings Kq which are proportional to the squares of odd integer numbers, and 
stationary densities on the first bifurcating branch are stable at least for small enough 
Kq. We discuss some of the results of previous authors in the light of our exact 
calculations. The Appendices discuss technical matters related to the bulk of the article. 



2. Probability density and moment-generating function 

In this section we consider general aspects of the model (|I|). First of all, we use 
the moment approach considered in to derive a nonlinear Fokker-Planck equation 
(NLFPE) for the one-oscillator probability density in the limit of infinitely many 
oscillators. The moment approach exploits the symmetry of the dynamical problem we 
are interested in, so is a good starting point to deal with the Kuramoto model which has 
rotational symmetry (an extension of this method to deal with tops has been considered 
in [|14|)- Secondly, we characterize phase and frequency synchronization generalizing to 
T 7^ the concept of order function introduced by Daido [§] for oscillator synchronization 
at T = 0. Lastly, we show that, at high frequencies, the probability density can be 
decomposed into m components rotating steadily at the frequencies of the peaks of a 
given multimodal natural frequency distribution. Each component probability density 
solves a NLFPE with zero natural frequency (in the rotating frame) and a modified 
coupling constant. This latter results follows immediately from the method introduced 



in Ref. |T7] for the usual Kuramoto model with a sinusoidal coupling function. This 
result can be combined with the exact solutions obtained in later Sections to yield 
analytic expressions of synchronized phases in models with frequency disorder and white 
noise forcing in the high-frequency limit. 
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2.1. Nonlinear Fokker- Planck equation 

To start with let us define, 

1 ^ 



^r = T^E (exp(^fc0,))^-, (3) 
i=i 

where the brackets denote average with respect to the external noise and the overbar 
average with respect to the random oscillator frequency. This set of moments is invariant 
under the local symmetry (pi <f)i+2TT which is the symmetry of the dynamical equations 
|1]. The equation of motion for the moments reads 

f)ffm oo 

— f = -K,tk Y: anH^+nH\ - k'THj: + tkH^K (4) 

^'^ n=-oo 

It is easy to derive from Eq.(^) the following differential equation 

— = -—\v(x t)g\ +T— - -^-^ (5) 

for the generating function 

g(x,y,t) = - Y: EexpHM^i/rW- (6) 

A:=-oom=0 

The velocity v{x,t) is defined by 

oo 

v{x,t) = -Ko E anH-nexp{inx) (7) 

n=— oo 

and it can be equivalently written as a convolution between the coupling and the 
generating function (see [Appendix A] ): 



v{x,t) = -Ko r f{x')g{x - x',0,t) dx'. (8) 
Notice that this drift velocity is related to Daido's order function H{x,t) (see 



Appendix B|) : 

v{x,t) = -KoH{x-net,t). (9) 

The case analyzed by Kuramoto corresponds to the force f{x) = sinx. Then ai = 
a*_i = —i/2 and the rest of components are zero. This yields v{x,t) = Korsm{6 — x), 
where = rexp(i6'). Let us assume that there exists a probability density p{x,LJ,t) 
such that 

1 ^ 

T7E-^(<^i(^)'^i) = / / J^ix,uj) p{x,uj,t)p{uj)dujdx, (10) 

iV J —IT J —OO 

as ^ OO. In particular, we assume that 



1 ^ 

p[x,uj,t)p{uj) = —Y^i^ - (Pjit)) - ^j) , (11) 
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tends to a smooth function in the hmit as — oo. We can relate the probabihty 
density to the moment-generating function as follows: 

1 oo CO m 1 A'' 

fe=— OO m=0 j=l 
1 Af 1 oo 1 ^ 

N 2n N^. ^^^^ ' ' 

j=l «;=— oo j=l 

e^"^ p{x,u;,t)p{uj)duj. (12) 



Then the probability density can be written as follows 

1 

2^ 



p{x,uj,t)p{uj) = ^ I g{x,iy,t)e dy (13) 

oo 



and 



r e-'^'^dy = -up{x,u;,t)p{u;). (14) 

By using these expressions in (^, we derive the usual nonlinear Fokker-Planck equation 
(NLFPE) for p{x,uj,t) wherever p{uj) ^ 0: 



dt ~'~ dx 



{u; + v)p=T^. (15) 



Equations (|T0]) and (P) provide us with the following formula for the drift velocity v in 
terms of p: 

rTT roc 

v{x,t) = —Kq / / f{x')p{x — x',uj,t)p{uj)dx'duj. (16) 



The probability density p should moreover satisfy a normalization condition J^^ pdx = 1 
and an initial condition p{x,u!,0) = pq{x,uj). A derivation of these problems by path 
integral methods can be found in ||TB[|. (The path integral derivation is applicable to each 
Fourier mode of the coupling function; see page 676 of []TB|)- See |TP|, ^ for rigorous 
proofs of ([TUp, ( P^ - (|I^) in different models. 

2.2. Phase and frequency distributions 

To characterize oscillator synchronization, it is convenient to define the phase and 
frequency probability densities. The phase probability density, P(</), t), is the probability 
of finding an oscillator with angle in (0, + d(f)) independently of its frequency. It is 
given by 



7r /-oo 



P{(f),t) = / / 6{x — (f)) p{x,uj,t) p{uj) dxdui. (17) 

J — IT J — oo 

Similarly, we may define the frequency of a given oscillator as the average (when it 
exists) 

1 r d(,y. 
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for sufficiently large r > 0. When N oo, (|T]) and the ergodic theorem imply that the 
previous time average is equal to 

1 /"^ 
CO -\ — / v{x,t)dt, 



T Jo 

where u = ujj. Then we may define the frequency density as 

/TT /"OO / 1 rr \ 

/ 6in-uj-- v{x,t')dt']p{x,n,t)p{uj)dxdn (18) 
-IT J-oo V T Jo J 

that can be simplified in the important cases of: stationary (i) and rotating-wave 
(ii) probability densities. In case (i), v = v{x) coincides with its time average, and 
is proportional to Daido's order function with Qe = 0. In case (ii), p{x,Q,t) = 
p{x — f2et. A), with A = Q — Qf, and, according d), v{x,t) = —KqH{x — Qgi). Then 
the change of variable ip = x — Q^t reduces probability density and order function to 
time-independent functions so that (|18[) becomes 



/TT roo 
/ 6[uj -Qe- A + KoH{ij)]p{ij,A)p{A)dijdA, (19) 



where p(A) = p{Qe + A). Daido considered only cases (i) and (ii) at T = with an 
order function having a single minimum Hmin (resp. maximum Hmax) oX if) = ipi (resp. 
= 4^2) HI- Then the probability density is 

A 
Ji'o^ 

+- — [i-x(^i,v>2)Wj . (20) 



p(^,A)=5(^-/f-i (A^) X{^.,^^m 

C{A) 
A-KoH{ij) 



where 



CjA) = \ - , (21) 

ip = (j) — flg^t, with A as before and X(ipi,i)2)i'^) equals 1 if ipi < ip < iIj2 and otherwise. 
On the interval —it < ip < it the order function H{ip) is stationary as before. C{A) is 
the frequency at which oscillators with angle outside {ipiji'2) rotate. Inside {ipi,ip2) "we 
have KqH{%P) = A and dip/dt = 0. Inserting (pO|) in (|T^) and (§) we find a functional 
equation for the order function, which can then be solved exactly or approximately 
Daido's expressions for the angle and frequency densities are found by inserting (pOD in 
(p!7|) and (p!9[) (minor notational changes have been made) |^ 

P{tP) = Kop{KoH{tP)) H\iP) + fl - X(^i,V2)(V')" 
^ C{A)p{A) 
A-K^HiiP) 

Pico) = Sico - n,) / piA)dA + i ^ . (23) 



/•°° C{A)p{A) 

J-oo A- K Hilb) ^(-^°°'^0Hrmn) (^) X(KoHma.,Oo) (A) rfA, (22) 
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2.3. High-frequency limit 

Further general considerations can be made for multimodal natural frequency 
distributions in the limit of high frequencies [l^. Let us assume that p{uj) has m 
maxima located at ujqQi, I = 1, . . . ,m, where Uq oo, and p{uj) duj tends to the limit 
distribution 

m 

r(fi)dfi = (24) 
1=1 

with 2^ a; = 1, and Vt = — , 

1=1 ^0 

independently of the shape of pioj) as ujq — ^ oo. pioj) doj and V{yL)dVL may be used 
interchangeably when calculating any moment of the probability density [including of 
course the all- important velocity function (0), which is related to Daido's order function 
as said before]. Thus any frequency distribution is equivalent to a discrete multimodal 
distribution in the high-frequency limit. The discrete symmetric bimodal distribution 
considered in [0, |^ corresponds to m = 2, 17^ = (—1)', ai = \, I = 1,2. By following 
the procedure explained in |jT^ , we can show that the oscillator probability density splits 



into m components, each contributing a wave rotating with frequency VLiuoq to the order 
function: 



p{x,u,t) = Y.pf\x-niuJot,t)+0{uj^^). (25) 
1=1 



The densities pf'^ obey the following Fokker-Planck equation 



m r^-/<><g^K' -5'W°'(/J'.<W'|-0. (26) 

where j3 = x — QiUot and each fZ-^pf^df] = 1. This equation corresponds to the NLFPE 
(p!5|) - (1^) in the moving variable (3, with a coupling constant K^ai instead of Kq and 
p{u!) = 6{uj). We shall see below that its solution evolves to a stationary state as the 
time elapses provided /(x) is odd. Thus the probability density (to leading order in 
1/uJo) is the sum of m components obeying the stationary solution of (^) with variables 
X — QiUJot. The overall velocity (related to the order function) is a superposition of m 
waves each traveling with frequency Qiujq 

v{x, t) = -Ko f: ai r fix - niuot - f3')pf\f3')d(3' , (27) 
1=1 

where pf\x) is the stationary solution of (p6|). 

We shall now consider NLFPE without disorder, i.e., ojj = or p{uj) = 6(00). Our 
results will be applicable to models with a multimodal natural frequency distribution in 
the high-frequency limit. For a model without disorder, the moment-generating function 
is independent of y and it equals the probability density: g{x, y, t) = p{x, 0, t) = g{x, t). 
If detailed balance is obeyed (i.e. if / is an odd function) then the model is purely 
relaxational and consequently the formalism of statistical mechanics can be applied. 
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3. Stationary solutions of models with odd coupling 

We can find functional equations for the order function of stationary or rotating wave 
solutions at nonzero temperatures without making Daido's assumptions. In this section 
we shall assume detailed balance, which occurs when / is an odd function extended 
periodically outside (— vr, vr) and satisfying f{x + tt) = —f{x). (The expressions for the 
general case are somewhat more complicated). Let us define 

V{x,t) = r v{s,t)ds, (28) 

J — TT 

W{x,LJ,t) = {n + x)uj + V{x,t). (29) 

The property /(x + vr) = —f{x) implies that the drift ([T6|) satisfies f (x + vr, t) = — f (x, t), 
and V{x + vr, t) = V{0, t) - V{x, t). 



3.1. Stationary solutions 



Let us restrict ourselves to the case of stationary solutions; rotating wave solutions may 
be reduced to this case after moving to a rotating frame. The stationary solutions should 
have the form 



p{x,u) 



Z e T 



J 



exp 



W{x,uj) -W{s,uj) 
f 



ds. 



(30) 



where Z and J are functions of u, independent of x. We now impose the condition that 
p be a 27r-periodic function of x, use the symmetry properties of the drift, and find the 
probability flux J as a function of Z: 



J 
f 



2 sinh 



T J 



ZZ^e — -T^dx 



(31) 



Z can be found from the normalization condition J^^ pdx = 1. The functional equation 
for the drift v{x) (or, equivalently, the order function) is obtained by inserting (|5DD, (pi|) 
and the formula for Z in (^B[). The result is 



ViXj 







fix 



X 



Z uj) 



2sinh 



27rtt; 



I e T — ds 

X — w(7:z-) 1 \ p{uj)dujdx\ 



(32) 



Ziuj) 



- w-(...) f2sinh(%^)/^',e- 



W{a,uj) 

T ds 



fl^e-^ds 



— 1 > dx. 



The functional equation (|3^) for v{x) may in general have several solutions, depending 
on p{uj) and the value of the parameters and T. Notice that no assumption on the 
shape of the order function has been made in order to derive (0). Thus we have a 
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general equation for v{x) valid for any temperature thereby extending Daido's theory. 
If f{x) is not odd, the same procedure yields a more cumbersome equation which we 
omit. An interesting case corresponds to the case without disorder, p{uj) = S{uj), for 
which J = 0, p{x, 0, t) = g{x, y, t) = g{x, t), and 



Vo(^) 

e T 



(33) 



The subscript zero just reminds us that we are considering stationary solutions. In this 
case there is a general Liapunov functional related to the free energy [^. In Sections 
^, 1^ and ^, we show that simpler quadratic Liapunov functionals may exist for specific 
forms of f{x). 



3.2. Liapunov functional 

Let us define the relative entropy 



g{x,t) 
dfi 

~dt ~ 



g{x,t) In 



V{x.t) — iJ.(t) 

e T 



dx. 



, , dV{x,t) , 
at 



(34) 

(35) 
(36) 



V"(x, t) is given by (pSf ) with the drift calculated by means of the exact probability 
density g{x,t). (p^) implies that 



d 

g{x,t) — \ng{x,t) dx = 0. 



Direct calculation shows that 
drj 



dx \g{x,t)^ 



dx < 0. 



-T / g{x,t) 

J — TV 

The inequality ?/ln?/>y — 1, Vy>0, can be used to show that 

/•TT 

^ > 1 — / ^(a;, t) dx. 



(37) 



(3^ 



Then the relative entropy is bounded below if e^*^^'*-* G L^(— vr, vr) and the function ji{t) 
is bounded below. The first condition thanks to the Trudinger-Moser theorem (see p^ ) 



is always fulfilled if V{x,t) is in the Sobolev space VFo^'^(— vr, vr) and this property is 



verified for V{x,t) defined in 
us write Eq. (^6]) as follows 
d_ 



). To show that the function /i(t) is bounded below let 



gV dx 



- V 

I —TT 



vg 



T 



dg_ 

dx 



V^dx 
9t 



dx = — v{x, t) J{x, t) dx, 



(39) 
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after using the nonlinear Fokker-Planck equation and integration by parts. J = vg 
—Tdg/dx is the probabihty flux. Then (|39|) can be equivalently written as 



= / [g{x, t) V{x, t) - g{x, 0) V{x, 0)] dx 

t rir 

/ v{x, s) J(x, s) dx ds. (40) 

Let us first note that, since J^^gdx = 1 and \^ is a uniformly bounded function, the 
two first terms in the right hand side of (^) are bounded. If we can prove that 

t r-K 

v{x, s) J{x, s) dx ds 

is uniformly bounded, then will remain bounded for all time and the demonstration 
that ?7 is a Liapunov functional will be over. By using the symmetry of f{x) and the 
definitions of J and v, it is easy to show first that 



r J{x,t)dx = 0, Vt. 

J — IT 



Since J is a continuous function with respect to x, the mean value theorem applied to 
the above equahty implies the existence of at G [— vr, tt] such that J{at,t) = 0, where 
the subscript makes reference to the fixed time t and at is chosen so that the function 
J has a definite sign over {at,x). Then 



/ / v{x,s) J{x,s)dx ds <m.a.'Xx^t\v{x,t)\ / signJ(a;, t) / J{x,s)ds 

Jo Jat ' Jao [ Jo 

where t = is the inverse function of x = a^. Now we integrate 



dx, 



dg_ ^ 

dt dx 



over {ao,x) x (O,^^), and obtain 



/ [9{y,t)-g{y,0)]dy=-T f f ^^^^^ dyds = -T f J{x,s)ds, 

Jat Jo Jat Oy Jo 

where we have used again Fubini's theorem to exchange the order of integration in the 
right hand side. The left hand side of this equation is bounded by 2, so that we obtain 
after taking the limit as t ^ cxd, 

ft I fx 2 

lim / J{x,s)ds < - lim / [g{y,t) - g{y,0)]dy < -. 



Combining this bound with the uniform one for v{x,t) and taking into account that 

E 



rt rn rt rx 

/ / v{x, s) J{x, s) dx ds = / v{x, s) J{x, s) dx ds, 

Jo J -IT jrt JO Jai 
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where A is a finite set of indices because of J G has a finite number of zeros in [— vr, vr], 
we deduce that the third term in ( ^0]) is also uniformly bounded and, as consequence, 
fi{t) is bounded from below. 

Our Liapunov functional may be used to show that g{x,t) tends to a stationary 
solution as time elapses and also to discuss the global stability of these solutions. Let 
Qooix, t) be the limiting probability density as t ^ oo. Equating the entropy production 
iW^) to zero, we obtain g^o = Pit) expV{x,t)/T. Inserting this into (p^Sf) , we find 

dp (3 dV{x,t) _ 
dt T dt 

which may be rewritten as 

Tdp _ dV{x,t) 
~p~dt ~ di ■ 

The left side of this expression is independent of x whereas the second side is not, so 
that both sides are zero. Then V is time independent and /? is a constant, which proves 
goo to be of the form 



3.3. Equilibrium states 

The previous considerations suggest that for models with odd-coupling functions and no 
frequencies [i.e. f{—x) = — /(x) and p{uj) = S{uj)], a thermodynamic formulation will 
suffice to identify the stationary states as well as the possible existence of thermodynamic 
singularities (i.e., bifurcations). 

In order to demonstrate this assertion, we shall start by defining an appropriate 
energy function. 



^ = ^T.Ei<l)^-<l^=) (41) 

i<j 

where E{x) is a two-pair interaction energy function defined by E{x) = J!^^f{s)ds and 
E{x) = E{x + 2n). The 27r-periodicity of E{x) is crucial to derive the results of this 
section. Note that this newly introduced function is related to the potential V{x,t) in 
Eq. (HD by 

_Y_^:^lJ. = I I E{x — x')p{x' ,uj,t)p{uj)dx'duj 

-tin J — oo J — TT 

/•TT 

E{—n — x')p{x' , u, t)p{uj)dx'duj. (42) 



— oo ^ — TT 



which has the physical meaning of an averaged energy. The second term is an additive 
constant which fixes the origin of the energy scale. Consequently the Liapunov functional 
defined in Eq. (|3^) is a generalized free energy for the system. Now we want to show 
that the equilibrium state of such a system yields potential solutions identical to (|33|). 
Let us compute the partition function, 

/"TT /"TT ^ 

Z= / .../ e-^'^Wdcl^, 



i=l 
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N 



exp |^-^E^(<^^ (43) 

We have Oq = because / is odd. Then we can write E{x) = I]^-oo bn^^""^ (plus 
an unessential constant term), where 6o = 0, 6„ = —i{an/n) {n ^ 0). Consequently we 
can rewrite the partition function for an oscillator system without frequency disorder, 
p{uj) = S{uj), as 

We have neglected the contribution of the terms {i = j) in the exponential (which is 
equivalent to redefine the origin of energies so that E{0) = 0). Let us now insert delta 
functions in the previous integrals and use the identity, 

= '•-'^p {"» (''» - ^ E ^"*) I '"'»''^»' 

where the new set of Lagrange multipliers A„ has been introduced. Inserting these 
representations of the delta function in (^) and permuting the orders of the integration 
between the (p's and the hn, A„, we reduce the final expression to a single site problem, 

/oo 
exp[NA{X,h)] H dXndhn, (45) 
-oo „ 



where 



A{X, h) = - J2 ^nK ^ Yl ^nKh^n 



n=—oo ^ n=— oo 

/ oo \ 



+ log / exp E A^e*"-^ (46) 

-'-^ \n=-oo / 

The dominant contribution to the partition function is determined by the saddle 
point equations, 

(47) 



dXn dhn 
which yield, 

hn= < exp^incj)) >, (48) 
A„ = - ^/i-„(&„, + &-„) = -f3Koh_^ Re(6„) . (49) 
Here the average < B{(f)) > of an arbitrary function B{(f)) is defined as follows 

" /-'.exp(E~_.A„e-"*)# ■ <^°' 

The next steps are quite standard. Inserting the expression (^91) for A„ into (H), we 
obtain a closed set of equations for the parameters hn- The reader will easily convince 
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himself that the parameters /i„ are nothing less that the moments defined in Eq. (|^) 
(with uji = 0). Finally we find that the effective Hamiltonian in the exponent of ( pOD is 
exactly given by the potential function Vo(x) of eq.(^). Then the equilibrium density 
function g{x) should satisfy 

/ \ u j.\ exp(/3V"(a;)) 

g{x) =< 6{x - 0) >= ^-^ — -^^^ , (51) 

J^^exp{pV{x))dx 



with 



V{x) = —Kq ^ {Rebn)h_nexp{inx)dx. (52) 



This coincides with the definition (^) [V{x) = y(a:)+constant] if all the 6„ are real i.e. 
if f{x) is an odd function. From (^S[) and ( ^Uj ) it follows that the function v{x) is a 
solution of the following functional equation 

vix) = -Ko < fix - 0) >= -Ko j^^^pvw ^0 • (53) 

The system of equations (^) and ( |5B| ) may have one, many or no solutions. In case of 
multiple solutions, the free energy and the dynamics should be used to fully characterize 
the equilibrium state. 

4. A family of models 

When p{uj) = 6{uj), Eq.(|D becomes 

^ = -^Kx,t),J+T— . (54) 
The probability density g{x,t) obeys the natural normalization condition 



g{x, t) dx = 1 . (55) 

Let us consider the following special cases of localized coupling functions /: 
(a) The hard-needles model: In this case f = 6' and (^) and (^ yield 

dg K ^ ( ^3\j^j'^^3 
dt dx^ dx' dx^ 

Note that in this case there is an energy of the model given by ii^ = Z]i<j (^(^i — <i>j) 

which is trivially zero except when two phases coincide. If we imagine the phases 

as needles located in the centre of the unit circle then the crossing of needles 

costs infinite energy. This is the reason for its name: In this case the needles 

are impenetrable and interact like the hard spheres in the theory of liquids. We 

shall see later that synchronization for this model is not possible at any temperature 

T > 0: any initial configuration g{x^{f) evolves towards incoherence [g = (27r)~^ 

which means that all angles have the same probability to occur] for large enough 
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time. This result is anticipated by considerations of linear stability alone. In fact, 
linearizing (^) about {2tt)'^ we obtain the equation 

d'g 27rA . _ 

for g[x, t) = {2ti)~^ + g{x) e'^*. This equation has 27r-periodic solutions of zero mean 
if —A = T+Kq/ {27r). Thus incoherence is linearly stable for any T > 0. In addition, 
as we will see in the next section, there exists a Liapunov functional which shows 
that g{x,t) should evolve in 7r,7r) towards x independent functions which, by 
normalization, is equal to (27r)^^. 

(b) The stick-needles model: Now / = —6' and we obtain the following evolution 
equation for g: 

— = -ir — f —)+ T— (57) 
dt ^ dx^ dx^ dx"^ ' 

The only difference with the previous case is the change of sign in the velocity field. 

This is enough to dramatically alter the behaviour of the model. Now the energy 

is given hj E = — J2i<j ^{(pi ~ 'Pj)- Obviously the thermodynamics is badly defined 

because the ground state energy is — cx). The needles want to stick to each other and 

there is only one relevant configuration which dominates the partition sum. Linear 

stability considerations indicate that incoherence is stable only if T > Kq/{2it). In 

general this model is not mathematically well-posed unless T — Kog{x,t) > 0. 

(c) The Burgers model: If f = 6, the velocity becomes v{x,t) = —K^g and the 
dynamical equation for g is 

dq do ^(9^0 , , 

^ = 2Kogjf + T-^. 58 
ot ox OX'^ 

The model corresponds to the Burgers equation (BE). Note that the two coupling 

functions / = ±5 are both equivalent (to go from the + to the - case or viceversa it 

is enough to make the transformation x —x.). The dynamics of the model 

corresponds to the Burger equation with a supplementary periodic boundary 

condition for the g (which plays the role of the velocity field in the BE equation). 

Performing the transformation x ^ x' = —x/2Kq we obtain the Burgers equation 

dg ^ ^dg _ ^d'^g ^^^^ 
dt dx dx"^ ' 

with viscosity u = jj^. Physically this model corresponds to a system of needles 







which tend to move together in the same direction when they meet. The incoherent 
solution is always a stationary solution of the BE equation. Straightforward but 
tedious calculation shows that the incoherent solution is the only stationary solution 
satisfying periodicity and normalization conditions. Notice that the coupling 
function for the Burgers model is not odd and therefore the general functional 
( p4D of Section 3 may no longer be a Liapunov functional. However, the functional 
G{t) = ^ J^^ g{x,t)'^dx satisfies G'{t) = — T J^^{dg{x,t)/dxydx, and is therefore 
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a Liapunov functional for T > 0, but not at zero temperature. Thus any 
initial configuration evolves towards incoherence for T > 0, while in principle 
synchronization is possible only at zero temperature. These results agree with 
considerations of linear stability for the incoherent solution: the probability density 
g = l/(27r) is linearly stable for T > and neutrally stable for T = 0. Finally, we 
shall mention that this model, contrarily to what happens in models (a) and (b), 
lacks a thermodynamic formulation because its dynamics violates detailed balance. 

(d) Daido coupling (the extended Burgers model): If we choose f{x) = sign(x) 
(periodically extended outside the interval [— 7r,7r]) as coupling function, the 
expression (|^) for the velocity becomes 

/n 
sign(0 gix- ^,t) di 
-IT 

= Ko\r g{x - e, t)di - r g{x - ^, t)d^ 

J-TT JO 

= -Ko r [gix -^,t)-gix-^ + tt, t)] d^. (60) 
Jo 



Then we have 

dv 



-2Ko[gix,t)-g{x + 7r,t)]. (61) 



The moment-generating function satisfies (0) and (|60D which form a nonlocal 
equation for g{x,t). Synchronization already appears at non-zero temperature 
(T > 0) for this simple model which is purely relaxational (when all the oscillator 
frequencies are equal to zero). The corresponding Hamiltonian is given by 

n = ^Y.^od{\<P,-<l),\,2n). (62) 

i<j 

In what follows we are going to analyze the dynamical behaviour of his family of models. 

5. The porous Medium models 
5.1. The hard-needles model 

The dynamical equation for the hard-needles model is 

dt ° dx y dx J dx"^ 
which is considered together with periodic boundary conditions 

, ^ , . , dg dg 
g{-TT,t) = g{n,t) and ^(~^' = ^(^' ^) (64) 

and the initial condition 

g{x,0) = go{x). (65) 
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5.1.1. Well-possed problem. This equation is a particular case of a the general quasi- 
linear equation in divergence form studied by O. A. Ladyzenskaja, V. A. Solonnikov and 
N. N. Ural'ceva in [2C], Theorem 6.1 and by S. N. Kruzhkov in |2^. As a consequence 
of the results in and we can deduce that the problem is well-posed if 



T > 0, 



which is always fulfilled for (7 > 0. 



5.1.2. Decay in time estimates. In order to obtain some decay estimates in time for 
the solutions of solution (|63D towards its equilibrium state, we must obtain a Green 
function to the linear heat equation 



dt dx^ 



(66) 



with periodic boundary conditions (|6^). Let T{x,t) be the heat kernel, i.e., the 
fundamental solution of the heat equation in IR, defined by 



r(x,t) = {A-nTty 



2 g 4Tt 



Let us recall the definition of the Theta function 28 , which is related to one of Jacobi's 



elliptic functions P9| , 

00 

e(x,t) = r(x + 2n7r,t), 

n=— 00 

for t > 0, and O = for t < 0, and define the function ^ by 



00 QY 



(67) 



(68) 



The functions G and \& satisfy 

e(x + 27r, t) = e(x, t) and ^(x + 27r, t) = ^(x, t). (69) 

Then, the function Q{x — z,t — s) verifies: (i) it is a solution of the linear heat 
equation ( pB] ) except in (s, z); (ii) it satisfies (p^); 6(x — z, t — s) — T{x — z,t — s) is 
solution of (1^). As a consequence, 9(x — z,t — s) is a Green function for the linear 
heat equation with periodic boundary conditions. Another interesting property of the 
function is its positivity. Also, note that the G function can be written as follows 



e(x,t) 



1 

2^ 



- V e~" ^* cos nx. 

" n=l 



In view of the preceding properties, we can write the solution of the hard-needles 
model in the following integral form 



g{x, t) = e{x - z, t)go{z) dz 

J — TT 

dQ{x — z,t — s) ( dg' 



+ 



TT JO 



dz 



^dz I 



(70) 



Exact solutions and dynamics of globally coupled phase oscillators 



18 



We will use for —7 > — ^ t > and x > the following estimate (see 



9^r(x,t) 



< 



c{Tt)- 



\x 



1-27+fc ' 



A; = 0,l,2,..., 



(71) 



where c = (47r)~^/^(— 476""'^ + 2e^^Y^~'^^''^'^ and we have chosen —7 such that —1/2 < 
—7 < and close enough to —1/2. 

Using these estimates, the generating function and its derivative with respect to 
the space variable can be estimated from (ffOD for small initial data go G L\^, where L^t, 
denotes the space of 27r-periodic functions belonging to Lp{—h, tt), see |^ for definition 
and main properties. This procedure was introduced in by G. H. Cottet and J. 
Soler to study the decay properties of the Navier-Stokes equations with weak initial 
data (singular filament measure). Let t G [0,r], r e R. Firstly, due to the positivity of 
0, go and of the solution, from equation (|63| ) we deduce 



go{x) dx. 



(72) 



To obtain the estimates of the solution let us note that since 6 converges 
uniformly we have 

1 

2^ 



0(x — z, t)go{z) dz 



go{z) dz 



00 



n=l 



cosn(x — z) go{z) dz. 



In the case of the uniform bound for g{x,t), we estimate (|70D as follows 



I-n9o{x)dx 



+ 



dx 



271 

t-s] 



< 



0(-,t) 



1 

2^ 



5'0 Li 

Ztt 



Ztt 



dg_^ 

dx 



ds. 



(73) 



Also, we have 



9{;t) 



I-n9oix)dx 



271 









< 








^■2 



\9o\\lI 



+ 



* dQ dg 
"a^("'^"")"^2JI^(-'")ll^§°JI^(-'")II^L^^' 



(74) 



and 



dg ^ 
dx 



X 



dg_^ 

dx 



< 



ZTT 



^1 

dx 



9o n + 

Ztt 



de^ 

dx 



LSI 



\9[Sr)\\Lsz 



Ztv 



d'g, 



dx"^ 



ds. 



rP/2 



(75) 



To close the circle of these estimates requires to prove that |^ G L^-k ■ This can be 
obtained directly from the equation (|^) combined with an estimate of |j G L2i^ deduced 
also from (^). We bound the previous inequalities using (|7TD . 
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Let Mp(t) be, with 1 < p < oo, such that 



{Tty 



l + C(Tt)l-27 

{Tty 

1 + C(Tt)l-27 



t) 



2tt 



go{x)dx 



< Mp{t). 



Z7V 



< Mp{t), 



Ztt 



Thus llff tends to as t ^ +oo if we choose —7 < —1/3. 

Set M = mlx{Mp{t), 1 < p < 00, t e [0,r]}. Then, we have for (|73|), and 
(fTsP the quadratic equation 

M <C{go) + C{T)M^, (76) 

where C{go) is a constant depending on the norm of go in Ll^. Inequahty ( |75D also 
imphes at least local existence for any initial data or global existence of the solution for 
small initial data in Ll^. For small initial data in Ll^ and t G [0, r], Vr G R, ( |TBD yields 



9{-,t) - — / go{x)dx 
/vr 



dg_^ 

dx 



t) 



< M{Tt)-^ 
< M{Tty (l + C{Ttf~'^^ 



C{Tt) 



1-27 



(77) 
(78) 



There exists other methods to study the existence properties but they do not give the 
above sharp estimates ([77|)-(|78D, see and p7[ . 



5.1.3. Asymptotic behaviour as t 00. Let us obtain that the system is simplified 
asymptotically as t ^ 00. In fact, we will prove that the solution of equation (^) 
converges in L^^ as t ^ 00 towards a limit function by using a Liapunov functional 
associated to the system. The functional 

G{t) = l r g{x,tydx, (79) 



(80) 



is positive and it satisfies 



G'it) = - 1 V + Kog) (^^^ dx < 0. 



Thus it is a Liapunov functional and g{x,t) should evolve towards x independent 
functions which, by normalization, equal (27r)~^. Since we have developed in the general 
case the study of the properties of the Liapunov functional, we refer to Section ^ for the 
detailed arguments of the above assertions. 



5.1.4. Convergence towards the solution of the Porous Medium Equation asT ^ 0. We 
want to solve equation ( |OBD from another point of view, and study the behaviour of its 
solutions as T ^ 0: we propose an argument of splitting in time (with step At) between 
the heat equation and the Porous Medium Equation (PME). Let us first consider the 
periodic porous Medium problem 



dg 



PM 



dt 



Kr 



d_ 

' dx 



9 



PM 



dg 



PM' 



dx 



^1) 
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,-PAf/ 



r -'[-7r,t)=g'^''\7r,t) 
g''^\x,0) = g,{x). 



and 



dg 



PM 



dx 



i-n,t) 



dg 



PM 



dx 



(82) 
(83) 



The initial and initial-boundary value problems associated to this equation were 
first studied around 1950 by Zel'dovich, Kompaneets and Barenblatt (see and p4| ). 



These authors explicitly obtained a fundamental (self-similar) solution of (0) in IR: 

|2 



C 



\x\ 



12ti 



(84) 



where [s]+ = max{s,0} and C is an arbitrary positive constant. This solution has 
compact spatial support for every fixed positive time. There is a free-boundary or 
propagation front, given by 



t = {12Cy 



\x\ 



which separates regions where g^^^ > from those where g^'^^ = 0. (|^) presents 
corners (jump discontinuities in its first derivatives) on the free boundary. In the 1-D 
case, the study of these problems for the PME started with the results of Oleinik and 
coworkers in ||3^. See the surveys by A. S. Kalashnikov 



or by J. L. Vazquez fS^] for 
more information. However, to our knowledge, the periodic boundary problem has not 
been studied so far. 

It is known, see that for all initial data go G L^^ there exists a unique 

weak solution g^^'^ G C([0, oo); Ll^), which is not classical in general. This is illustrated 
nicely by the self-similar solution (p^, which is 27r-periodic. In our periodic boundary 
context, this selfsimilar solution makes sense, at least locally for < t < 7r7(12C)2, 
before the propagation front arrives to the boundary. This example illustrates that a 
solution (corresponding to arbitrary initial data) may have jump discontinuities in its 
first derivatives on the propagation front or free boundary separating the regions where 
the solution is positive from those where g^^^ = 0. Moreover, disturbances propagate 
with finite speed. In general, it is clear that g^'^^ > for t > 0, provided go > 0. 
Solutions of (PT|)-(|S1) have the following behaviour as t — oo: 



lim 

t— >oo 



t) 



I-^9o{x)dx 



2n 



0, VpG[l,oo]. (85) 
], by following the 



Lp{-tt,tt) 

This result can be proved, according to C. M. Dafermos's ideas 
basic steps listed below: 

1) Prove that the semigroup U corresponding to the problem ([HT|)-(P3|) is a continuous 
contraction semigroup on L^(— 7r,7r) which verifies the maximum principle 



||f^(^)5'o||LP(-7r,7r) < ||5'o||LP(-7r,7r) 



and 



go < go a.e. in [-7r,7r], implies Uit)go < U{t)gQ. 
2) The orbit 7 ((^o) = Ut>of/(t)(7o is relatively compact in L^(—7r, tt) for (^o G L°^(— 7r,7r). 
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3) The cu-limit set uj{uq) is non-empty and compact in L^. 

4) Use an appropriate Liapunov functional such as 

1/(0 = ess sup^g[_^^^]Oa;). ^ G ^^(-vr, vr), 

and the contractive property of the associated semigroup to prove that (i) V is a. 
constant W on uj{uo), and (ii) the cj-hmit set consists of constants. 

5) The comparison principle given in 1) and the fact that the average of the solution 
is preserved allows us to identify W with 



1 

2^ 



go{x) dx. 



and to prove (|^) for p = 1. 

The result follows for p > 1 by using the dominated convergence theorem. 

A similar proof based on Dafermos's theory was introduced by N. D. Alikakos and 
R. Rostamian in in the case of boundary condition of type 



dx 



(x) = 0, for X = —71, 71. 



We refer to their paper for the details of the above scheme of proof. 

It is also possible to estimate the rate of decay towards the equilibrium for the 
solutions. In fact, if go G L°°, there exist positive constants a and k such that (see EDI) 



PM/ ^^ I-n9o{x)dx 



g'^'''{;t) 



2n 



L°°(-n,n) 



We now come back to the complete hard-needles problem for nonzero temperature. 
The splitting method consists of solving in [0, At] the PME (T = 0) and diffuse it by the 
heat kernel (T > 0) to obtain gi = g^{x,At) = Q{At,x) *x U{At)gQ, where g^ means 
the splitting solution and is 

e(x, At) fix) = r e(x - y, At)f{y) dy. 

The method is well defined and follows the same scheme in every interval [nAt, {n—l)At]. 
The iterative map (?„ U{At)gn-i ^ Q *x U{At)gn-i is a composition of contractive 
maps in L^. The diffusivity of B implies g > 0, for t > 0, and the nonexistence of 
a propagation front, i.e. of a free boundary. Also the splitting argument shows the 
convergence in of the solutions of the hard-needles equation to the PME as T — 0. 

5.2. The stick-needles model 

The dynamical equation for the stick-needles model is 

— = -K — ( —\ + T— (86) 
dt ^ dx \ dx I dx^ 
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An analysis similar to the hard-needles model may be done in the present case. 
However, there is a hurdle to be resolved previously: the model is not well-possed unless 
(see 0) 

T 

Moreover this condition implies the non convergence to a well-possed problem as T ^ 
because the generating function must be non negative. 

6. The Burgers model 

dg dg d'^g 

- = 2K,g- + T— (87) 

The existence and uniqueness properties of this equation are well understood, see 
^ for the main results and references. In the Hopf-Cole transformation is used 



to analyze the Burgers equation with periodic boundary conditions. 

To deduce some some decay estimates for small initial data in L^^, we can repeat 
the same arguments as in the study of the hard-needles model by using the integral 
formulation 

g{x,t) = / &{x - z,t)go{z) dz 



1 dQ{x - z,t ~ s) 2/ 



This procedure leads, for 1 < p < oo, to 

\\9{-,t)-7r / 9o{x)dx\\LP < ct •^II^oIIli • (89) 

In the case of the Burgers equation it is also possible, as in the Porous Medium 
case, to deduce that 

G{t) = \ fj{x,tfdx 
is a Liapunov functional only if T > and its derivative satisfies 



(90) 



Then, repeating the same ideas as in Section ^, we can prove that g{x,t) should evolve 
towards (27r)~^ if T > 0, after normalization. At zero temperature G{t) is an integral 
of motion and the previous argument cannot be used. Figure^ shows the time evolution 
of G{t) for different temperatures. For T = it remains constant while for T > it 
decays towards the expected value. 

On the other hand, the spatially periodic solution of the Burgers equation converges 
towards an entropy solution of the Hopf equation 
dg" fjdg" 

-w = ''''' ^' ^''^ 
g"{-T^,t) = (?^(7r,t), (92) 

(7^(0, x) = go{x), (93) 
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Figure 1. Time evolution of the Liapunov function for the foUowing values of 
temperature: a) T=0 (dotted line), b) T=0.05 (dashed line), and c) T=0.1 (solid 
line). 



as T ^ 0. This follows by using same arguments as in the hard-needles model. In 
P^l , K. L. OuYoung proved the following result based upon ideas developed by M. G. 
Crandrall for the Cauchy problem in R: the semigroup S : go S{t)go 
contracting semigroup in Ll^ and we have 

^^St)|| 



g^ is a 



1^ 



< 



\9o\\L° 



a.e., 



TV{g''{;t))<TV{go), 

f \g''{t2,x)-g''{ti,x)\ dx<cTV{go)\t2-ti\, 

Jtr 

for ti, t2 > 0. The uniqueness of solution for ( pT]) is based in the fundamental work of 
S. N. Kruzhkov, 

Then, the splitting argument used in the hard-needles model give the convergence 
in Ll^ for the solution of the Burgers equation to the entropy solution of the Hopf 
equation. 

It is well-known that in general (^Tj), ( p2D and ( p3D do not have continuous solutions 
for all time because shock waves appear after finite time and tend to dominate the 
solution as t ^ oo. The quantities J^^g^dx are preserved in time until a shock is 
formed. After a shock wave appears, only the mass J^^g dx is conserved along the 
time evolution. Figure ^ is an illustration of the phenomenon. E. Hopf in discussed 



Exact solutions and dynamics of globally coupled phase oscillators 



24 




Figure 2. Seven snapshots of the generating function g taken at different times. We 
have chosen as initial state a Gaussian distribution. Notice the appearance of the shock 
wave. 



the asymptotic behaviour for large time of solutions of (|9l|)- (|9^) . Later, P. D. Lax 
and T. P. Liu and M. Pierre |^ completed and generalized this study for an arbitrary 
conservation law. A first result in this direction is given by the following inequality 

271 



\g^{x,t) - M{x,t)\ < c- 
where M{x, t) is the mean function 



t 



(94) 



M{x,t) = — 
Ztx 



{x, t) dx. 



If go is in L^, then the mean function is time independent and holds. 

Let us assume that initially go{x) has a single maximum and a single minimum at 
each period. Then g^ evolves towards a precise asymptotic shape so that its profile 
for large time looks like a 27r-periodic sawtooth profile with a single shock per period. 
Between shocks g^ decreases linearly. The asymptotic behaviour profile result in the 
case of a Cauchy problem was specified by T. P. Liu and M. Pierre in By adapting 



the arguments in we can prove for our periodic problem 



lim t^p' \\g"i-,t) - uj{- 



t~*co 



Z7V 



0, 



Exact solutions and dynamics of globally coupled phase oscillators 



25 



for an initial data Qq G L^^ and 1 < p < oo. Here 

1 X 1 / Kq \ 

to (x,t) = — — — — = -— X 1 — n < X < n (95) 

^ ' ^ 2tt 2Kot 2Kot \ tt J ^ ^ 

which is a sawtooth profile moving towards the right at speed Kq/tt. The strength of the 
shock max((yf) — min(5f) = 7r/{KQt) converges to zero as the time goes to infinity. Thus 
the shock weakens, it vanishes and we are left with the incoherent solution (27r)~^. In 
summary, incoherence is the globally asymptotically stable solution for the periodic 
Burgers equation at any temperature, including zero. In the last case, interesting 
transient patterns including shock profiles are possible before every memory of the initial 
conditions is washed out. 

7. The extended Burgers model 

Instead of working with the nonlocal model given by Equations (|^) and (|60D for g{x, t) 
we can try to find local equations for f (x, t) and g{x,t). Let us define 

p{x,t) = 2Ko[g{x,t) - g{x + Tc,t)], 

a{x, t) = 2Ko [g{x, t) + g{x + vr, t)]. (96) 
These functions obey the equations (l60|) and 



dp _ d{av) ^ ^d'^p ^g^^ 
dt dx dx"^ ' 

dt dx^ dx-' dx"^ 



Now we can obtain an equation for v{x,t) by inserting (plD , i.e. dv/dx = —p in (|9 
and integrating the result with respect to x: 
dv d'^v 

--a.-T-^^m. (99) 

Equations ( |60D and (|96D together with the periodicity of g{x,t) imply that 

cr(x + vr, t) = cr(x, t), v{x + TT,t) = —v{x,t), (100) 

i.e. the 27r-periodic functions cr{-,t) and v{-,t) are periodic and antiperiodic in the x 
variable, respectively. Moreover the normalization condition for g and the vr-periodicity 
of a imply 

a{x,t) dx = 2Ko. (101) 





The conditions ( p.00| ) may be used to show that J{t) = in (|99D, so that v obeys 
dv r^d'^v 

g^-..-T^ = 0. (102) 

To see this, substitute x + vr as the argument of v and a in (p9D, and use ( |100| ) to obtain 
J = — J . Once cr and v are found, we obtain the generating function from the equation 

,(,.,).J_U.,)-M^1. (103) 
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A similar development to the study of the existence properties of solutions given 
for the Porous Medium and Burgers models can be done for v{x, t) in ( p,02| ). In fact, we 
can write 

6(x — z, t)vQ{z) dz 

-TT 

+ r Q{x - z,t - s){av){z,s) dz ds (104) 

J-n Jo 

coupling it with 

/TT 
6(x — z, t)ao{z) dz 
-TT 

rTT rt QQU - Z,t - S) f dv\ , , , , 

+ L i ^ifc — ' j <^»^' 

The same type of estimates as in the previous sections together with an iterative method 
and a fix point theorem lead to a cubic equation, instead of a quadratic one as in ([76|) . 
This allows to ensure the existence and bounds in L^, at least for small initial data in 
L^, of both v{x, t) and a{x, t), for positive temperature T. However, the Theta function 
does not appear with a space derivative in (|104|) as in the Porous Medium and Burgers 
models which implies the non-convergence in time towards its respective mean values. 

For this extended Burgers model the asymptotic in time behaviour is described by 
the Liapunov functional analyzed in Section ^ Let us now find the stationary solutions 
of these equations. From (^) and (|102|) , we obtain 
^d'^v 

Ta + — = L + Mx. 

a{x) and v{x) are 27r-periodic functions obeying ( |100| ) and L and M are constants. 
27r-periodicity of a and v imply that M = in the second of these expressions. Thus 

Ta + ^ = L. (107) 

Now we have L > because a > 0. Combining ( |101| ) and ( |107| ), we find the following 
equation for L: 

2TKq+ r — dx = TTL. (108) 
JO 2 

Inserting (|107|) into (|106|) we obtain 
This may be integrated once yielding 

E 




2T^ 8T^ ' (110) 

The energy E, < E < L^/2 is calculated so that the period of v{x) be 2tt. It is 
convenient to rewrite the previous equation in terms of 
/.X v(x) 
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Then ( |1 10| ) becomes 

1 / dw\^ 

2 \di) 
Suppose that we choose w{fS) ■ 
Then w{^) is 

r^i^) dw 



w 



E 
1? 



= £. 



(112) 



2£ 



0, w'(0) > as initial conditions for a given trajectory. 

(113) 



until V reaches the turning point wq{£), where 

Wo{£) = V2 sJl-VT^^ , (114) 

at ^ = P{S)/A. Notice that we will obtain a one-parameter family of stationary solutions 
because + c), with w{^) as in ( |1 13] ) and c G R+, is also an admissible stationary 
solution. Because of symmetry, an oscillation period is first completed at four times this 
value, i.e., when 

rMS) dw 4:y/2K{m) 



P{S) 



1 + VI - 2£ 



m 



1 - vr^2^ 



;ii5) 



(116) 



1 + VI -2^ 

Here K{m) is the complete elliptic integral of the first kind with parameter m given 
by ( |116| ); see ||2^, page 590. P{S) is an increasing function from P(0) = 2tt to 
P(l/2) = +00. Branches of admissible stationary solutions satisfy ( |100D , in particular 
they fulfil v{x + vr) = —v{x). Then their energies £ G (0, 1/2) should be such that 

27rVI 



nP{S) 



T 



for positive odd integers n = 2p — 1, p 
inadmissible solutions such that f (x+vr) = 
Eq. ( |108| ) for L may be written as 

rwoie) w'^dw 

n 



(117) 

= 1,2,.... Notice that even n's would yield 
: +v{x). By changing to the w and ^ variables, 

vrVI 2Ko 



28 



T 



[n odd) or, equivalently, as 



Vsyi + Vl-2^ [K{m) - E{m)\ 



ttVI 



2Kn 



118 



(119) 



T VI 

Here E{m) is the complete elliptic integral of the second kind with parameter m given 

by (inl). 

Given i^'o > fixed. Equations ( p,15|) , ( [L17| ) and ( |119|) yield the values of £ and L 
corresponding to a branch of stationary solutions with a given odd value of n. A given 
stationary solution {n,S,L) will have n maxima and n minima. We can find L from 
( |117| ) and eliminate it in ( |119| ) with the result 

An^ K{m) 



E(m) 



VI - 2SK{m) 

1 + vr^2f" 



T 



(120) 
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This latter equation relates the energy S to the coupling parameter Kq. As Kq increases 
so does the number of possible stationary branches. Notice that the stationary branch 
with index n exists for values of the coupling constant larger than Ko = n^T7r/2. At 
this value, S = (therefore, m = 0). Clearly the number of possible stationary branches 
is then half the integer part of ^J2Kq/ (ttT) (only odd n are admissible), which increases 
with Ko. For a given stationary solution v{x), the stationary probability density may 
be reconstructed from (|103|) as 

]_ iL_f__dv) 
T ~ 2T ~ dx\ ' 



'1211 



A stationary periodic moment-generating function is given by (|110|) , (|1 1 1|) , (|113|) , 
( |117|) , (|118|) , ( |120|) and ( |121| ). Their explicit form may be obtained by using the 
definitions of the Jacobi elliptic functions; see ^9], pages 569 and ss. From (|111|) , 



( |TT3|) , ([ml) , ([T20D and (|), we find 

AnTy/mKim) 

v[x) = 
H{x) -- 



■ SUM . 



TT 



n 



ma 



u 



E{m) 
2nK{m)x 



\/l-2£ K(m) 



snu , 



TT 



The probability density may be reconstructed from ( |1 1 7| ) , ( |120| ) - ( |1 24| ) : 



K{m) 
27r 



1 



msn M 



m CUM dnw 



Eim) - 



(122) 
(123) 

(124) 
(125) 



Notice that these synchronized solutions are defined up to a constant phase shift 
X ^ X + c, as we said before. Such solutions exist for any positive temperature T > 
if Kq > n'^TTT/2. Therefore, it is possible to find synchronized solutions of the extended 
Burgers model for any positive temperature. The stationary drift velocity ( p.22|) has 
n = 2p — 1 [p = 1,2, . . .) maxima and n minima on the interval — tt < x < ti. Between 



successive extrema, v{x) vanishes once. According to (|121|) , the extrema of g{x) on the 
interval (— tt, vr] are reached at the zeros of v{x). A little algebra shows that 



dx ~ AKqT"^ 



and thus 



dx"^ 



v=0 



AKqT^ 



w dw 



dw 




Since dw/d^ = ±v2£^ wherever w = 0, we have 



'^^'^ d^ 



v=0 



( dw^ 
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Therefore maxima of g{x) are reached at points where v = and v'{x) < 0. These zeros 
are u = 2lK{m), where / is an odd number running from —n + 2 to n. The corresponding 
values of x are = li^/n. Thus an oscillator population described by (|122|) is split 
in n subpopulations with angles close to x = a;;^„, / = — n + 2, . . . ,n. The frequency 
density is found by means of (^8|). In dimensionless form, p„(fi) = \/L Pn{VLQ,), which 
yields 



Here the sum is over all solutions ^ G (— vr, vr] of the equation w{C,) = fi. The maxima of 
Pn{^) yield the likeliest frequencies of rotation for the oscillators. See [|, ^ for additional 
physical interpretation and numerical simulations. 

Once a stationary solution (with index n), g{x) [equivalently we may specify 
cr(x) and v{x)], has been found, it is interesting to discuss its linear stability. Let 

t) = cr(x)+e'^*o"(x), v{x, t) = v{x)+e^''v{x). Then a and v satisfy (|10CI|) , Jq adx = 0, 
and solve: 

d'^v 

-T—-[a{x)v + v{x)a] = -Xv, (127) 

^d'^a d'^[v{x)v] ^„ (12S) 
dx"^ dx^ 

Solving these equations in the general case seems difficult. These equations are easy to 
solve for the incoherent solution v{x) = 0, a{x) = 2Kq/'k. In this case (|127[ ) and ( |128| ) 
uncouple and we explicitly find the eigenvalues 

X = ^-(2p-lfT, (129) 

TT 

A = ~{2pfT, (130) 



(where p = 1, 2, . . .) corresponding to (|127|) and ( |128|) , respectively. Eq. ( |129| ) shows that 
incoherence is linearly unstable for Kq > T'k/2. At Kq = {2p — l)^T7r/2 the stationary 
solutions (|125|) (with n = 1, 3, . . .) bifurcate from the incoherent solution. This can also 
be checked as follows: linearize (|109|) about v = and calculate the value of Kq for 
which there is a 27r-periodic solution of the resulting equation which satisfies ( p.00| ). It 
is Kf) = Ttt/2 corresponding to = 0, n = 1. This solution is linearly stable (at least 
for values of Kq in an appropriately small half- neighbourhood of T'k/2) because of the 
principle of exchange of stabilities between bifurcation branches. In figure |^ we show 
the different bifurcation branches which occur at different integer values oi n = 2p — 1. 

An alternative calculation of the first branch (index n = 1) of stationary solutions 
bifurcating from incoherence is using Monte Carlo methods. They are quite powerful 
since acceptance rate can be tuned at will in such a way that the approach to 
the stationary state may be faster. Furthermore, the intrinsic dynamics in Monte 
Carlo methods is discrete in time, unlike the continuous-time dynamics of the original 
problem. Thus the rounding errors due to time discretization are absent in Monte Carlo 
calculations of stationary solutions. We use the Monte Carlo method and the Glauber 
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Figure 3. Bifurcation diagram £ versus 2Ko/{ttT) for stationary synchronized 

solutions branching off from incoherence 5 = at the square of the odd integer 
numbers. 



algorithm with a random sequential updating of the phases for the model Hamiltonian 
(|6^) of Section |[(d). Local phases are randomly changed from (pi to (pi + aS, where 
a takes values with equal probability on the interval [0, 1] and S is the typical size of 
the move. The proposed change is accepted with probability [1 + exp(/3A7Y)]^^ where 
ATi is the change in the Hamiltonian. A convenient synchronization order parameter 
r is defined through the global mag netization M = re*° = (1/A^) E^i e*'^^ The order 
parameter can be calculated from Monte Carlo simulations of the Hamiltonian (^) with 
E[x) = \x\ and Kq = 1; see figure Notice that the curves in this figure tend to a 
curve intersecting the horizontal axis at the bifurcation point T = 2/7r as iV increases. 
This corresponds to the expected bifurcation result for the NLFPE (dot-dashed line in 
Figure | corresponding to = oo). Notice that finite-size corrections are of order A^~^ 



far from the bifurcation point and of order N~i near it ||T9 . 

The transition temperature corresponding to the first branch is easily obtained 
through standard finite-size scaling methods. Consider the kurtosis (or Binder 
parameter) for the synchronization parameter r, i? = |(3 — ^^2^2 ) ? where < . . . > is 
the standard configurational average [weighted with the usual Boltzmann-Gibbs factor, 
exp(— and Ti is given by (^)]. The curves for B are shown in figure ^ for different 
sizes. Note that these curves (specially for = 50, 100, 500, data for A^ = 1000 is more 
noisy) intersect at a common point characterizing the bifurcation temperature. 

It is interesting to compare the present results with those of Daido's order function 
theory 0. Except for a constant factor, —Kq^, Daido's order function is just v{x,t) 
in our notation. Notice that as Kq/T 00, ( |123D and ( [125| ) corresponding to the 
linearly stable solution with index n = 1 become H{x) = — sign(a;) and g{x) = 6{x — tt), 
respectively. This coincides exactly with Daido's solution at T = lll|, P]. An interesting 
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Figure 4. Synchronization parameter r as a function of T for different sizes 
N = 50, 100, 500, 1000, oo. 




0.40 0.50 0.60 0.70 0.80 0.90 



Figure 5. Binder parameter B as a function of T for different sizes (in the low T 
region, from bottom to top) N = 50, 100, 500, 1000. The curves become steeper for 
larger values of A''. The theoretical value of the bifurcation temperature T = 2 /it is 
marked in the figure. 
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aspect of our exact construction of stationary solutions is that they can shed some 
hght on scahng near bifurcation points as T ^ [jl^, ^. In fact, ( |120| ) imphes that 
stationary solution branches issue forth from incoherence at Ko = k„ = n^7rT/2 (n odd). 
Furthermore ( |120| ) shows that the nth synchronized stationary solution branch satisfies 



£ oc {Kq/T — n n /2) for Kq/T close to its bifurcation value Kq = Kn- Then ( |122| ) - ( |125| ) 
show that g{x) = 0{\Kq — nd^^'^) if T > stands for the value of the bifurcation 
parameter Kq/T at the bifurcation point Clearly as T ^ 0+, a quasicontinuum of 



= 0{\J Kq/T) stationary branches has bifurcated from incoherence for a small fixed 
Kq. In these circumstances, the derivation of a one-mode amplitude equation as in Ref. 
||TT| does not describe correctly the situation. See for a derivation of an amplitude 
partial differential equation describing a quasicontinuum of Hopf bifurcations; similar 
techniques could be used in the present case. 

The results at zero temperature may be obtained in another form. The model with 
Daido coupling is described by (0), (|102|) and ( |103|) . Setting T = in these equations, 
we obtain an integro-differential equation 



dt Jo dx \ dx 



which can be consider as a modification of the Porous Medium equation with memory. 
This equation can be also written in a local form in the following way: 

I dv da 1 ^ 

a 



V dt' dt 2 dx"^ 



Then v{x, t) obeys 

d (ldv\ _l 9V 
Wt V^Stj ~ 2 'do^ 
d^u 1 S^e^" 



(131) 

dt^ 2 ' (132) 

if \v\ = e". This is an amplitude equation with infinitely many terms (e.g., expand the 
exponential in a Taylor series), and g{x,t) = 0{Kq) in agreement with Daido's results: 
the critical exponent is (3 = 1 ii g = 0{\Kq — Kc\^). Nevertheless, (|131|) or (|132|) are 
very different from Crawford's amplitude equations . 



Clearly the stationary solutions of (|131|) compatible with the symmetry requirement 
v{x + vr,t) = —v{x,t) are Vq{x) = q sign(x), where g is a constant. In particular, 
incoherence corresponds to g = 0. Cq^x) obeys av = 0. Thus ctq is zero except at the 
points where f = 0, namely x = 0,n (mod. 2tt). If ( p.01| ) holds and a{x + n) = a{x), we 
have 

ao{x) =2Ko[6{x) +6{x-tt)], (133) 

on the interval x G (— vr, vr]. Inserting dvo/dx = —2q [S{x) — 6{x — 'Jt)] and (|133|) in (|103|) , 
we obtain 

goix) = 6ix-n), (134) 



Exact solutions and dynamics of globally coupled phase oscillators 



33 



provided q = Kq. This density function is the same as that obtained above by taking 
the hmit Kq/T ^ cxd of the stable stationary density with index n = 1. Dai do found the 
same results by using his order function theory p. Notice that the stationary solution 
( |134| ) is a member of the one-parameter family of stationary solutions, gQ{x + c), c = 
constant. 

Linearization of (^) and ( |102| ) (with T = 0) about Vo{x) = Kq sign(a;) and (|133| ), 

v{x, t) = Vq{x) + V{x, t), a{x, t) = ao{x) + S(a;, t) {V <^ Vq, T, <^ ao), yields 

dV 

— = ao{x)V + Vo{x)J:, (135) 
9S dvoix)V 

which is equivalent to 

W{x,t) = vo{x)V{x,t). (138) 

If we prove that W{x,t) evolves towards zero as e~^*^, M > 0, then this implies that 
the density ( 134|) is linearly stable. But analyzing Equation ( 137 ) as a Heat equation 



with respect to t, with boundary conditions W ^ as \t\ —>■ oo, we have the Green 
function of the linear part 

T{x,t) = (47rii'o"^|a;|)~2e , x > 0. 

Using the periodicity and symmetry properties of v and a, the solution can be written 
in convolution form with respect to t with initial data KQV{0,t). This allows to assure 
that W{x,t) evolves towards zero as e~*^*^, with M = Kq/Att. 



8. Concluding remarks 

In this paper, we have considered the Kuramoto model for synchronization of phase 
oscillators with global periodic coupling functions f{x) and subject to external random 
forces. For general couplings and natural frequency distributions, we have derived the 
one-oscillator nonlinear Fokker-Planck equation and given an approximate formula for 
its solution in the limit of high natural frequencies. Provided the frequency distribution 
has m peaks in the high-frequency limit, this formula indicates that the one-oscillator 
probability density splits into m components. Each component corresponds to the 
solution of a NLFPE with zero natural frequency, on a frame rotating with fixed angular 
velocity. Which results do we know for such a reduced zero-frequency NLFPE? 

In Section ^ we have shown that a Liapunov functional of the free energy type exists 
for a class of odd coupling functions (zero natural frequency). Then the probability 
density evolves towards a stable stationary function. Which stationary function this is, 
depends on f{x). We have analyzed a family of singular coupling functions for which 
the NLFPE reduces to partial differential equations or systems thereof, such as the 
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porous media equation, the Burgers equation or systems of Burgers equations. For 
these equations, we have examined the behavior of their solutions in the hmit of large 
times for different temperatures; see Section ^ for details. The most interesting coupling 
function, f{x) = sign x (periodically extended outside [— 7r,7r]), was already studied by 
Daido in the case of zero temperature, T = (see Section |^.2 for a rephrasing of his 
assumptions and results). This model is capable of synchronization at any temperature. 
We have found exact formulas for stationary synchronized probability densities which 
bifurcate supercritically from incoherence. Although we have not been able to determine 
their linear stability, we know that the first bifurcating branch of solutions is stable 
at least for coupling parameter near the bifurcation point. Let us assume that this 
bifurcating branch is always stable for coupling parameter above the bifurcation value. 
We can combine our results to obtain the stable probability density for a model with 
a multimodal natural frequency distribution in the high-frequency limit. We find that, 
except for a constant shift in x, the stable probability density is 

m 

P-EP^l^-^^^ot), (139) 

where 9{x) is the Heaviside unit step function. In these equations, g{x; K^ai) is the 
solution ( |125| ) corresponding to setting n = 1 and KqUi instead of Kq in (|120|) . The 
overall velocity function (proportional to Daido's order function) is a superposition of 
rotating waves correponding to the contributions of the synchronized components 



v{x, t) ~ l^aiO [ Koai - — ^ 

1=1 ^ 2 y 71 

X su pf(.n.,)(.-n,.„t)j ^^^jj 
Here rrii is the value of m obtained by substituting n = 1 and KqUi instead of Kq in 

O- 

Appendix A. A convolution identity for the velocity 

To obtain the convolution identity (|^), we can write the velocity v{x,t) as follows 

oo 

v(x,t) = -Ko «ni^V"" 



'5 

n=— OO 
oo 



- ^0 E r e-^-^'fix')dx') Hly 

- Ko— / fix') E e<^-^'^H\ dx' 

27r J-TT 



which, using the expression for g{x,t) 
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gives the announced relationship (||). 



Appendix B. Order function and oscillator drift velocity 

In our notation, Dai do's order function is defined as 

oo 

H{x,t) = - hkZke~'^\ 

k=—oo 

where Daido's coupling function is related to ours by h{x) = — /(— x) {hk = —a-k), and 

1 ^ /.TT / 1 \ 

Zk{t) = ^ ^^e^M0.W-ne*] ^ J^^^^k. I Ys{<P,{t) - Q^t - x) \ dx 
= r e'''^g{x + Qet,0,t)dx. 

J — n 

Here we have used that A^~^ ^f=i ^i'Pji^) " ^) tends to g{x, 0, t) in the limit of infinitely 
many oscillators. Thus Zk{t) is 2tt times the Fourier coefficient of g{x + flet,0,t) with 
index —k. Inserting this in H{x,t) above, we find 

H{x,t) = - J\{x' + net,0,t) i £ hke"'^'''-Adx' 

\k=—oo I 

= - r g{x' + nj, 0, t) h{x' - x) dx' = T g{x' + Vt^t, 0, t) f{x - x') dx' 

J —TT J — 77 

/TT 
g{x + figt — x', 0, t) f{x') dx', 
"TT 

which, together with (H), implies 

H{x, t) = --^ v{x + ilet, t). 

Then (^) follows, as said in the Introduction. 

Acknowledgments The authors want to express their gratitude to Rafael Ortega 
and Juan L. Vazquez for useful discussions and for pointing them some references. We 
acknowledge partial support by the DOES (Spain) Projects PB98-0142-C04-01 (LLB) 
and PB98-1281 (JS), FOM (The Netherlands) contract FOM-67596 and DOES (Spain) 
contract PB97-0971 (FR) and TMR (European Union) contract ERB FMBX-CT97-0157 
(LLB & JS). 



References 

[1] A.T. Winfree, Geometry of Biological Time, Springer- Verlag, New York (1990). 

[2] Y. Kuramoto, Self-entrainment of a population of coupled nonlinear oscillators, in International 

Symposium of Mathematical Problems in Theoretical Physics, H. Araki ed.. Lecture Notes in 

Physics, Vol.39 (Springer, New York), (1975) 420-422. 
[3] S. Shinomoto and Y. Kuramoto, Phase transitions in active rotator systems, Prog. Theor. Phys. 

75, (1986) 1105-1110. 

[4] Y. Kuramoto and I. Nishikawa, Statistical macrodynamics of large dynamical systems. Case of a 
phase transition in oscillator communities, J. Stat. Phys. 49, (1987) 569-605. 



Exact solutions and dynamics of globally coupled phase oscillators 



36 



[5] Y. Kuramoto, Chemical Oscillations, Waves and Turbulence, Springer, Berlin (1984). 
[6] S. H. Strogatz, Norbert Wiener's brain waves, in Lect. N. Biomath. 100, edited by S. Levin, 
Springer, N. Y. 1994. 

[7] K. Wiesenfeld, P. Colet, and S.H. Strogatz, Synchronization transitions in a disordered Josephson 

series array, Phys. Rev. Lett. 76, (1996) 404-407. 
[8] H. Daido, Order function and macroscopic mutual entrainment in uniformly coupled limit-cycle 

oscillators, Prog. Theor. Phys. 88 (1992), 1213-1218. 
[9] H. Daido, A solvable model of coupled limit-cycle oscillators exhibiting partial perfect synchrony 

and novel frequency spectra, Physica D 69 (1993), 394-403. 
[10] H. Daido, Onset of cooperative entrainment in limit-cycle oscillators with uniform all-to-all 

interactions: bifurcation of the order function, Physica D 91 (1996), 24-66. 
[11] J. D. Crawford, Scaling and singularities in the entrainment of globally coupled oscillaiors, Phys. 

Rev. Lett. 74 (1995), 4341-4344. See also J. D. Crawford and K.T.R. Davies, Synchronization 

of globally- coupled phase oscillators: singularities and scaling for general couplings. Physica D 

125 (1999), 1-46. 

[12] L. L. Bonilla, C. Perez- Vicente, F. Ritort and J. Soler, Exactly solvable phase oscillator models 
with synchronization dynamics, Phys. Rev. Lett. 81 (1998), 3643-3646. 

[13] C. J. Perez- Vicente and F. Ritort, A moment-based approach to the dynamical solution of the 
Kuramoto model, J. Phys. A (Math. Gen.) 30 (1997), 8095-8103. 

[14] F. Ritort, Solvable dynamics in a system of interacting random tops, Phys. Rev. Lett. 80 (1998), 
6-9. 

[15] J. D. Crawford, Amplitude expansions for instabilities in populations of globally- coupled oscillators, 

J. Stat. Phys. 74 (1994), 1047-1084. 
[16] L. L. Bonilla, C. Perez- Vicente and R. Spigler, Time-periodic phases in populations of nonlinearly 

coupled oscillators with bimodal frequency distributions, Physica D 113 (1998), 79-97. 
[17] J. A. Acebron and L. L. Bonilla, Asymptotic description of transients and synchronized states of 

globally coupled oscillators, Physica D 114 (1998), 296-314. 
[18] L. L. Bonilla, Stable Probability Densities and Phase Transitions for Mean-Field Models in the 

Thermodynamic Limit, J. Stat. Phys. 46, (1987) 659-678. 
[19] D. A. Dawson, Critical dynamics and fluctuations for a mean-field model of cooperative behavior, 

J. Stat. Phys. 31 (1983), 29-85. 
[20] P. Dai Pra and F. den Hollander, McKean- Vlasov limit for interacting random processes in random 

media, J. Statist. Phys. 84 (1996), 735-772. 
[21] L. L. Bonilla, J. C. Neu and R. Spigler, Nonlinear stability of incoherence and, collective 

synchronization in a population of coupled oscillators, J. Stat. Phys. 67 (1992), 313-330. 
[22] L. L. Bonilla, C. J. Perez Vicente and R. Spigler, Time-periodic phases in populations of nonlinearly 

coupled oscillators with bimodal frequency distributions, Physica D 113 (1998), 79-97. 
[23] L. L. Bonilla, J. A. Carrillo and J. Soler, H-theorem for electrostatic or self-gravitating Vlasov- 

Poisson-Fokker-Planck systems, Phys. Lett. A 212 (1996), 55-59. 
[24] L. L. Bonilla, J. A. Carrillo and J. Soler, Asymptotic behavior of an initial-boundary value problem 

for the Vlasov-Poisson-Fokker-Planck system,, SIAM J. Appl. Math. 57 (1997), 1343-1372. 
[25] J. M. Kim and J. M. Kosterlitz, Growth in a restricted solid-on-solid model, Phys. Rev. Lett. 62 

(1989), 2289-2292. 

[26] O. A. Ladyzenskaja, V. A. Solonnikov, N. N. Ural'ceva, Linear and quasilinear equations of 
parabolic type. Translations of Mathematical Monographs 23, AMS, Providence, 1968. 

[27] S. N. Kruzhkov, Periodic solutions to non-linear equations, Differential'-nye Uravnenija, 6 (1970), 
731-740. 

[28] G. B. FoUand, Introduction to PDEs. Princeton U. P., Princeton, 1976. 
[29] M. Abramowitz and LA. Stegun, Handbook of Mathematical Functions. Dover, N. Y., 1965. 
[30] D. V. Widder, "The Heat Equation", Academic Press, New York-San Francisco-London, 1975. 
[31] O. Vejvoda, "Partial differential equations: time periodic solutions", Martinus Nijhoff Publishers, 



Exact solutions and dynamics of globally coupled phase oscillators 



37 



The Hague-Boston-London, 1982. 
[32] G. H. Cottet and J. Soler, Three-dimensional Navier-Stokes equations for singular filament data, 

J. Diff. Equations 74 (1988), 234-253. 
[33] Ya. Zel'dovich, A. S. Kompaneets, Theory of heat transfer with temperature dependent thermal 

conductivity, in Collection in honour of the 70th birthday of academician A. F. loffe, Izdvo. 

Akad. Nauk SSSR, Moscow (1950), 61-71. 
[34] G. I. Barenblatt, On some unsteady motions of a liquid gas in a porous medium, Prikl. Mat. Mekh. 

16 (1952) 67-78. 

[35] O. Olcinik, S.A. Kalashnikov, Y. L. Czhou, The Cauchy problem and boundary-value problems for 

equations of the type of unsteady filtration, Izv. Akad. Nauk SSSR, Ser. Mat. 22 (1958), 667-704. 
[36] A. S. Kalashnikov, Some problems of the qualitative theory of non-linear degenerate second-order 

parabolic equations, Russian Math. Surveys 42 (1987), 169-222. 
[37] J. L. Vazquez, An introduction to the mathematical theory of the porous medium equation, in 

"Shape Optimization and Free Boundaries", M. C. Delfour ed., Math, and Phys. Sciences, Series 

C, vol. 380, Kluver Acad. Publ., Dordrecht-Boston-Leiden (1992), 347-389. 
[38] L. A. Peletier, The Porous Medium Equation, in "Applications of Nonlinear Analysis in the Physical 

Sciences", H. Amann et al. eds., Pitman, London (1981), 229-241. 
[39] C. M. Dafermos, Asymptotic behavior of solutions of evolution equations, (M. G. Crandall, Ed.) 

Academic Press, New York, 1978. 
[40] N. D. Alikakos and R. Rostamian, Large Time Behavior of Solutions of Neumann Boundary Value 

Problem for the Porous Medium Equation, Indiana U. Math, J. 30 (1981), 749-785. 
[41] E. Godlewski, P. A. Raviart, Hyperbolic systems of conservation laws, SMAI 3/4, Paris, 1991. 
[42] A. Parker, On periodic solution of the Burgers equation: a unified approach, Proc. R. Soc. Lond. 

A 438 (1992), 113-132. 

[43] K. L. OuYoung, Periodic Solutions of Conservation Laws, J. Math. Anal. Appl. 75 (1980), 180- 
202. 

[44] S. N. Kruzhkov, First order quasilinear equations in several independent variables. Math. USSR-Sb. 
10 (1970), 217-243. 

[45] E. Hopf, The partial differential equation Ut -\- uux = l-iUxx, Comm. Pure Appl. Math. 3 (1950), 
201-230. 

[46] P. D. Lax, Hyperbolic systems of conservation laws H, Comm. Pure Appl. Math. 10 (1957), 537- 
566. 

[47] T. P. Liu and M. Pierre, Source-Solutions and Asymptotic Behavior in Conservation Laws, J. Diff. 

Equations 51 (1984), 419-441. 
[48] L. L. Bonilla and F. J. Higuera, The Onset and End of the Gunn Effect in Extrinsic Semiconductors, 

SIAM J. Appl. Math. 55 (1995), 1625-1649. 



